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ABSTRACT 

Field  trials  when  observations  are  correlated  with  those  in  neighboring 
or  nearby  plots  in  one  and  two  dimensions  are  analyzed  using  simultaneous 
autoregressive  models.  The  relationships  between  the  maximum  likelihood 
solutions  and  the  corresponding  well-known  Papadakis  estimators  are  clarified 
and  it  is  shown  that  the  maximum  likelihood  solutions  are,  for  specific  types 
of  designs,  easier  to  obtain  directly  then  by  iterating  on  the  Papadakis 
estimator  as  has  been  suggested.  A  simulation  study  compares  the  different 
models  and  methods  for  one  and  two  dimensions.  - 

Note i  In  spite  of  its  earlier  number,  MRC  Technical  Summary  Report  No.  2650 
is  a  sequel  to  this  report. 


AMS  (MOS)  Subject  Classifications:  62K10 

Key  Words:  Correlated  plot  yields,  Maximum  likelihood  estimator, 

Papadakis  estimators.  Simultaneous  autoregressive  field  model 
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SIGNIFICANCE  AND  EXPLANATION 

The  problem  of  analyzing  data  from  field  trials  when  observations  in 
neighboring  or  nearby  plots  are  correlated  is  an  important  one  which  has  been 
extensively  examined  over  the  years.  In  1937,  Papadakis  suggested  a  sensible 
but  mysterious  treatment  estimator  form  whose  role,  and  relationship  to  the 
corresponding  maximum  likelihood  estimator  has  been  debated  for  forty-odd 
years.  These  matters  are  here  clarified  for  certain  one-and  two-dimensional 
designs  under  certain  simultaneous  autoregressive  models.  Some  simulation 
results  are  also  given. 


The  responsibility  for  the  wording  and  views  expressed  in  this  descriptive 
summary  lies  with  MRC,  and  not  with  the  authors  of  this  report. 


ROLE  OF  THE  PAPADAKIS  ESTIMATOR  IN 
ONE-  AND  TWO-DIMENSIONAL  FIELD  TRIALS 

N.  R.  Draper  and  D.  Faraggi 
1.  INTRODUCTION 

The  analysis  of  field  trials  when  observations  in 
nearby  or  neighboring  plots  are  correlated  has  received 
extensive  study  over  a  number  of  years.  Early  work  was 
done  by  Papadakis  (1937) ;  his  suggested  method  of  model 
estimation  has  been  regarded  as  the  standard  procedure  and 
its  consequences  have  been  re-examined  by  a  number  of  authors 
including  Atkinson  (1969),  Bartlett  (1978),  Martin  (1982), 
Wilkinson,  Eckert,  Hancock  and  Mayo  (1983);  for  examples 
using  the  Papadakis  method,  see  Kempton  and  Howes  (1981) . 

It  has  been  pointed  out  by  Atkinson  (1969) ,  Bartlett  (1978) , 
Ripley  (1978,  discussion  to  Bartlett,  1978)  and  Martin  (1982) 
that  the  Papadakis  estimator  is  an  approximation  to  a  maximum 


likelihood  estimator.  The  contexts  of  their  remarks  differed 
somewhat.  Atkinson  (1969)  was  discussing  the  one-dimensional 
"plots  in  a  line"  case,  following  up  work  by  Williams  (1952) . 
Bartlett  suggested  iterating  the  Papadakis  procedure  in  one 
and  two  dimensions.  Ripley  (1978)  provided  a  matrix  expres¬ 
sion  for  the  maximum  likelihood  estimator  for  one  and  two 
dimensions,  and  showed  its  equivalence  to  generalized  least 
squares.  For  additional  discussion  see  Ripley  (1981, 
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pp.  88-97).  Martin  (1982)  also  used  a  matrix  formulation 
and  showed  that  Bartlett's  iterative  Papadakis  procedure 
converged  to  the  maximum  likelihood  estimator  in  certain 
circumstances . 

The  present  paper  considers  the  one-and  two-dimensional 
situations  in  matrix  form.  We  clarify  the  nature  of  the 
Papadakis  estimator  in  general  and  show  its  relationship  to 
the  maximum  likelihood  estimators  arising  from  two  types  of 
error  correlation  structure.  We  show  further  that  in  certain 
cases  maximum  likelihood  solutions  are  obtainable  sufficiently 
straightforwardly,  so  that  the  Papadakis  estimator  is  not 
actually  needed. 

In  general,  we  assume  all  one-dimensional  designs  to  be 
circular  in  the  sense  that  the  nth  plot  will  be  supposed  to 
be  to  the  left  of  the  first  plot,  as  would  be  the  case  if 
the  string  of  plots  formed  a  collar  on  a  mountain  peak.  A 
practical  approximation  to  such  an  arrangement  might  involve 
adding  additional  treatments  before  and  after  the  plot  string, 
setting  the  treatment  of  plot  n  in  a  "0  plot"  and  setting 
the  treatment  of  plot  1  in  an  " (n+1)  plot".  See,  for  example, 
Dyke  and  Shelley  (1976)  or  Draper  and  Guttman  (1980). 

Similarly,  two-dimensional  designs  are  assumed  to  be 
torus  designs,  wrapped  around  in  both  dimensions,  as  assumed 
by  Martin  (1982). 
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2.  ONE  DIMENSION,  ONE  SIDED  ERROR  STRUCTURE 

2.1  The  First  Order  Autoregressive  Model. 

We  suppose,  following  Williams  (1952)  and  Atkinson 
(1969),  that  an  experimental  design  consists  of  a  string  of 
adjacent  plots  used  to  examine  t  treatments.  Each  treat 
ment  will  occur  m  times.  If  each  treatment  occurs  c 
times  adjacent  to  every  other  treatment,  the  design  is 
called  a  Type  11(a)  design  (by  Williams,  1952).  If  c  *  2 
t  *  m  +  1.  Note  that  these  designs  sure  assured  to  be  circular. 

We  suppose  a  directional  association  exists  between 
errors  in  the  following  sense.  Let  there  be  n  response 
observations  yi'y2'***'yn  in  plot  sequence  from  left  to 
right.  The  assumed  model  is 

Y*  8  Og  ^  1  8  l,»».,n  s  *  l,...,t  (2.1.1) 

where  ag  denotes  the  effect  of  the  s*"*1  treatment. 
Furthermore 


xi  *  P*i^i  +  ei  (2.1.2) 

where  e  =  (e. ,e,, . . . ,e  ) *  ~  N(0,I c2)  . 
i  z  n 

A  first  order  autoregressive  relationship  be- 


tween  the  x's  is  assumed  instead  of  the  usual  indepen¬ 
dent  and  identically  normal  distribution  assumption.  The 


joint  distribution  of  xi'x2,,#*'xn  is 


f ^Xl/X2' *  *  * ,xn^  ~ 


n-1 

I  <«i> 

i*2  1 


(1  -  pM^exp^T  jx2  +  x*  +  (1  + 

- 2p  jx  (xixi+i>}  ] ; 


(2.1.3) 


see  Koopmans  (1942)  and  Box  and  Jenkins  (1976,  p.  276.) 

The  maximum  likelihood  estimator  for  p  can  be  found  in 
several  steps.  From  the  derivative  of  the  log  likelihood 
with  respect  to  p  we  obtain 


e: 


u-e2) 

a 


[n-l  n-i 

Jx  (xixi+i»  - ?  Jx  (xi>. 


n-l 


(2.1.4) 


which  is  the  maximum  likelihood  of  a2  in  terms  of  @  . 
Secondly,  from  the  derivative  of  the  log  likelihood  with 
respect  to  a2  and  a  substitution  for  o2  by  its  maxi¬ 
mum  likelihood  estimator  (2.1.4)  we  find 
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rn-1 

nlJlXiXi+l 


-  0 


X 


+  x2  + 
n 


d+a2) 


(2.1.5) 


which  can  be  solved  iteratively  for  $  .  To  find  the 
maximum  likelihood  estimator  of  the  effect  of  the  s*'*1 
treatment  ag  ,  substitute  x^  *  v^  -  ag  into  (2.1.3); 
the  Jacobian  is  the  identity  matrix.  This  gives 


-d2 


3L 


3a 


s 


u  +  02)  I 

Ci3s 


(yi-as) 


-  0 


L  <y.  - 


[itlj-s 


4ti:> 


(2.1.6) 


and  the  maximum  likelihood  estimator  of  a.  is 

s 


s 


a 

l+S 2 


l  <yA 

Ci±l]=s  1 


(2.1.7) 


where  the  w  acknowledges  Williams  (1952) . 

Note  that  y  y.  is  used  to  denote  the  sum  of  all 
[i]«s  1 

the  responses  from  the  plots  receiving  treatment  s  .  So 


O  «.*  ■.  n  '.  *.  *.'  »  •.*  *. v.  \  •.  •/•.  \  \  • 
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£  y.  is  the  sum  of  the  responses  from  the  plots 
ti±l]=s  1 

adjacent  to  the  plots  receiving  treatment  s  . 

2 . 2  The  Papadakis  Procedure 

Papadakis  (1937)  suggested  a  sensible  but  somewhat 
mysterious  estimation  procedure  described  by  Atkinson 
(1969) .  The  corrected  yield  of  the  i*"*1  plot  receiving 
treatment  s  is  defined  by 

y*  =  y  -  nf1  Y  y  (2.2.1) 

C  i  J=s 

The  concomitant  variable  6^^  is  defined  as 

«i  *  (yl-x  +  Vi+i)/2  •  (2.2,2) 

The  Papadakis  estimator  of  the  effect  of  treatment 


s  is  then 
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where 


n  * 

l 

i=l  1  1 


n 

l  «! 

i=l  1 


(2.2.4) 


and  p  acknowledges  Papadakis  (1937) . 

Atkinson  (1969)  showed  that  the  expectation  of  6 
when  using  a  large  sample  approximation  to  x^  and  taking 
the  expectation  of  the  ratio  to  be  the  ratio  of  expecta- 
tions  is  E(B)  *  2p/(l+p2).  Comparison  of  the  two  esti¬ 
mators  (2.1.7)  and  (2.2.3)  implies  the  Papadakis  estimator 
is  an  approximation  to  the  maximum  likelihood  estimator 

Ay  a 

oi  ,  with  a  . .  in  the  maximum  likelihood  estimator  corre- 
s  sti 

sponding  to  m  ^  \  y.  in  the  Papadakis  estimator, 

Ci]=s±l  1 

and  p/(l+p2)  in  the  maximum  likelihood  estimator  replaced 
by  6/2  in  the  Papadakis  estimator. 


i  * 


8 


2.3  The  Estimators  in  Matrix  Form 

Both  the  MLE  and  the  Papadakis  estimator  can  be 
written  more  conveniently  in  matrix  form.  Define  T  to 
be  the  design  matrix  of  size  t  by  n  whose  s^  row 
contains  1  in  the  i^  column  if  treatment  s  is 
applied  to  the  i^  plot  and  zeroes  otherwise.  For 
example ,  for  the  design  (1,2, 3,4) (2, 3, 1,4) (3, 1,4, 2),  the 
T  matrix  will  be 


T 

4x12 


"loo 
0  10 
0  0  1 
0  0  0 


0  0  0 
0  10 
0  0  1 
10  0 


10  0 
0  0  0 
0  0  1 
0  10 


1  o  (T 
0  0  1 
0  0  0 
0  10 


. (2.3.1) 


Then,  with  an  obvious  matrix  and  vector  notation. 


Y*  =  Y  -  m-1  T'  T  Y  =  (  I  -  m"1  T*  T  )  Y  (2.3.2) 
nxl  nxl  nxt  txn  nxl  nxn  nxt  txn  nxl 


where  Y*  =  (yj,yj, . . . ,yj) ' ,  y=  lY1'Y2' • • • 'Yn> '  so 
that  the  treatment  averages  are  m-1  TY  .  To  get  the  con¬ 
comitant  vector  6  we  define  the  neighbor-specification 
matrix  N  of  size  n  by  n  whose  l  row  contains  l  in 
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positions  j  for  which  plot  j  is  adjacent  to  plot  i 
and  zeroes  otherwise. 

In  our  example  it  will  be: 


N  = 
12x12 


0 

1 

0 

0 

0 

0 

0 

0 

0 

0 

0 

1 


1 

0 

1 

0 

0 

0 

0 

0 

0 

0 

0 

0 


0 

1 

0 

1 

0 

0 

0 

0 

0 

0 

0 

0 


0 

0 

1 

0 

1 

0 

0 

0 

0 

0 

0 

0 


0 

0 

0 

1 

0 

1 

0 

0 

0 

0 

0 

0 


0 

0 

0 

0 

1 

0 

1 

0 

0 

0 

0 

0 


0 

0 

0 

0 

0 

1 

0 

1 

0 

0 

0 

0 


0 

0 

0 

0 

0 

0 

1 

0 

1 

0 

0 

0 


0 

0 

0 

0 

0 

0 

0 

1 

0 

1 

0 

0 


0 

0 

0 

0 

0 

0 

0 

0 

1 

0 

1 

0 


0 

0 

0 

0 

0 

0 

0 

0 

0 

1 

0 

1 


1 

0 

0 

0 

0 

0 

0 

0 

0 

0 

1 

0 


(2.3.3) 


hen  6  =  i  NY* 
nxl  * 


and  the  vector  of  Papadakis  estimators  is 


aP  =  m-1  {TY  -  i  §TN(I  -  m-1T'T) Y) 
txl  1  ' 
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(2.3.4) 


where  aF 


ffi 


where  ar  =  *  *  *  ,at>' *  The  estimator  B  can  also 

be  expressed  in  matrix  notation  as  a  ratio  of  quadratic 
forms : 


2Y ' (I-m“1T'T)N (I-m‘1T'T) Y 
Y' (I-m“1T'T)N2 (I-m~1T'T) Y 


(2.3.5) 


The  maximum  likelihood  estimator  in  matrix  notation  is 


aW  =  m~ 1  { TY  -  -£^T  (TNY  "  TNT'aW)}  . 


This  implies 


aW  =  m_1{  (I  -  m"'l‘YTNT, )  ~A  (TY  -  YTNY)  ) 


”  ^vinxim  i  I 


=  (TV_1T' ) ~1TV'1Y  , 


(2.3.6) 


where 


a.W  (aW  £W  '  *WV  i 

a  ^1 ^ 9 


Y  =  5/d  +  P2)  , 


-1 

V  =  I  -  yn  . 


■-  •-  ■;  ■;  -.A,< 

■ ' '  s'a'  v 


Note  that  TT'  =  ml.  Equation  (2.3.6)  is  equivalent 
to  (23)  of  Wilkinson  et  al.  (1983) ,  as  pointed  out  by 
discussants  to  the  latter,  and  also  to  (5.34)  of 
Ripley  (1981) . 

It  is  well  known  that 

(I  -  A)"1  =  I  +  A  +  A2  +  A*  +  ...  (2.3.7) 

if  and  only  if  all  the  eigenvalues  of  the  matrix  A 
are  smaller  than  one  in  absolute  value.  We  can  thus 
apply  the  above  expansion  to  a  if  all  the  eigen¬ 
values  of  m_1YTNT'  are  less  than  one  in  absolute 
value.  The  matrix  TNT'  is  of  basic  importance  in 
these  sorts  of  analyses  because  its  (s,w)t*1  element 
is  the  number  of  times  treatment  s  appears  adjacent 
to  treatment  w  in  the  design.  We  shall  later  need 
and  define  matrices  of  the  form  TN^T'  to  express 
more  distant  spatial  relationship.  For  Type  11(a) 
designs  with  c  -  2, 
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m"1YTNT' 


a 

b 

b 


b 

a 

b 


>  b 
(m+1) 


b  •  •  •  •  b 

b  •  •  •  •  b 

a  •  •  •  •  b 

•  •  • 
•  • 
•  •  • 
•  • 

b  •  •  •  •  a 
x  (m+1) 


(2.3.8) 


A 

where  a  »  0,  b  =  2m  Y,  and  has  the  two  distinct  eigen¬ 
values  (Rao,  1973,  p.  67) 

1.  a  -  b  «  -2m  ^Y , 

2.  a  +  mb  =  2Y . 

Both  eigenvalues  here  are  less  then  one  in  absolute  value 
which  guarantees  the  convergence  of  the  corresponding  ex¬ 
pansion.  Expanding  gives 

oW  «  m“1{  (I  +  m-1YTNT'  +  m”2Y*TNT'TNT'  +  •  •  • )  (TY  -  YTNY)  } 
«  m“1{TY  -  YTNY  +  m~1YTNT'TY  -  m_1Y2TNT'T  +  •  •  •  } 

=  m"1{TY-YTN(I-m"1T,T)Y-m“1Y2TNT,TN(I-m”1T'T)  Y 

-  m"2Y3TNT'TNT,TN(I-m"1T,T)Y  -  •••}.  (2.3.9) 


Because  E(8/2)  ■  Y,  we  see  that  the  Papadakis  estimator 
picks  the  zero  and  first  order  terms  in  YN  from  the  auto- 


regressive  process,  while  ignoring  the  higher  order  terms. 
Note  also  the  repetitive  application  of  the  Papadakis 
correction  term  with  higher  and  higher  power  of  m-^YTNT ' 
in  the  full  expansion. 

Bartlett  (1978)  suggested  that  the  Papadakis  procedure 

should  be  an  iterative  one,  i.e.,  the  corrected  yields 

y?  88  y.  -  m”1  I  y.  should  be  changed  in  the  j  ^  iter- 
1  1  C i j=s  1 

ation  to  y?(^  =  yi  -  £P(j-l)  to  give 

ap(^  =  m_1{TY- j  EtN(Y  -  T'ap(^”1) )  }  (2.3.10) 

where  ap(0)  -  m-1T'Y  so  that  ap(1)  is  (2.3.4).  Martin 
(1982)  showed  that  indeed  this  procedure  converges  to  the 
maximum  likelihood  estimator  for  Type  XI (a)  designs,  for 
8  fixed.  Convergence  occurs  because  the  Papadakis  ex¬ 
pression  (2.3.4)  consists  of  a  truncated  portion  of  a 
convergent  expansion  (2.3.9)  of  the  maximum  likelihood 
estimator  (2.3.6). 

Some  useful  neighbor  balanced  designs  suggested  by 
Williams  (1952)  are  given  in  Appendix  1. 


2.4  An  Exact  Solution 


w 


In  fact,  for  Type  11(a)  designs,  iterative  use  of 


the  approximate  Papadakis  formula  is  unnecessary.  The 


(I-m  tTNT* )  matrix  has  a  form  similar  to  (2.3.8)  where. 


for  example,  for  type  11(a)  designs  with  c  =  2, 


a  =  1,  b  =  -2m-1Y  . 


(2.4.1) 


This  allows  it  to  be  inverted  as  a  patterned  matrix 


(Rao,  1973,  p.  67).  In  general,  the  inverse  of  such  a 


pattern  matrix  (see  Rao,  1973,  p.  67)  is 


(2.4.2) 


b'  b'  b* 


where 


a'  = 


a+ (n-2)b 


[a+ (n-l)b][a-b] 


(2.4.3) 


Ca+(n-l)b][a-b] 


«?? 


This  provides  the  exact  solution  to  2.3.6 


mCl-2y][l+2?/m]  {l +  m=2('n^l}?  Yrt*r  Hty-yTOY} 


(2.4.4) 


Substitution  of  this  in  Equation  (2.1.5)  provides  a 
single  equation  to  solve  for  p,  e.g.,  via  the  Newton- 
Raphson  method.  The  solution  is  then  used  in  (2.4.4). 
(The  existence  of  this  explicit  solution  method  seems  not 
to  have  been  noted  by  previous  authors . ) 


3.  ONE  DIMENSION,  TWO-SIDED  ERROR  STRUCTURE 


3.1  The  Simultaneous  Autoregressive  Model 


If  we  apply  the  one-sided  error  structure  approach 
of  Section  2  to  the  two  dimensional  case,  geometrical  im¬ 
balance  occurs.  This  can  be  rectified  by  a  two-sided  error 
structure  assumption,  as  we  shall  discuss  in  Section  4. 
Once  this  point  is  appreciated,  it  is  then  natural  to 
return  to  the  one-dimensional  case  and  consider  the  con¬ 
sequences  of  a  two-sided  error  structure  assumption.  We 
thus  assume  (following  Besag,  1974,  simultaneous  auto¬ 
regressive  model,  p.  201)  that 


yi  “  as  +  xi  i  *  8  ”  (3.1.1) 

where,  now. 


xi  ■  p(xi-i  +  xi+i)  +  €i  (3-1-2) 

with  e  ~  N (0 , la2 )  .  The  likelihood  function  is  (Besag, 
1974,  p.  201,  Equation  4.13  taking  y  =  0) , 
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(2Tro2)"n/2|B|exp{-(2o2)‘1(X'B'BX)}  ,  (3.1.3) 


where 


B  =  I  -  pN  , 


X  -  Y  -  T'a, 


(3.1.4) 


and  N,  T,  Y  and  a  are  defined  as  before.  The  maximum 
likelihood  estimator  for  a  is 


a  *  (TV_1T' ) -1TV_1Y, 


(3.1.5) 


where 


V-1  **  (I  -  pN) 2 


(3.1.6) 


Note  that,  since  TT'  »  ml,  (3.1.5)  can  be  rewritten  as 


and  since 


6t  =  m_1{[ I-m” *T ( 20N-0 2N2)T']-1TCl-( 20N-0  2N2 )  ]Y }  (3.1.7) 


N2  »  N(2) 


+  21  » 


(3.1.8) 


(2)  . 
where  N  is  the  "lag  two"  neighbor-specification  ma- 

th 

trix  whose  i  row  contains  1  in  positions  j  for  which 
plot  j  is  adjacent  but  one  to  plot  i  and  zeroes 
otherwise.  An  alternative  form  is 


3.2  The  Modified  Papadakis  Estimator 

The  expansion  of  the  inverse  matrix  [I-m  1T(2pN-p2N2) 
T']_1  in  (3.1.7)  converges  if  and  only  if  all  the  eigen¬ 
values  of  m~^T (2pN-p2N2 ) T'  are  less  then  one  in  absolute 
value.  For  example,  for  designs  of  Type  III  (Williams, 
1952)  which  are  obtained  by  imposing  simultaneously  two 
conditions  on  the  design: 

(i)  That  each  treatment  should  occur  equally 
often  adjacent  to  every  other  treatment. 

(ii)  That  each  treatment  should  occur  equally 
often  adjacent  but  one  to  every  other 


treatment. 
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For  example,  the  design 


(1,2, 3,4)  (2, 1,3, 4)  (1,3, 2, 4) 


(3.2.1) 


has  each  treatment  occurring  twice  adjacent  to  every  other 
treatment,  and  twice  adjacent  but  one  to  every  other 
treatment  assuming  the  design  is  circular,  that  is  the 
final  4  is  adjacent  to  the  initial  1.  The  matrix 
m"1T(2pN-p2N2)T'  is  of  form  (2.3.8)  with 


a  «  -2p2 ,  b  =  (4p-2p2 ) /m  . 


(3.2.2) 


Again,  such  a  matrix  has  two  distinct  eigenvalues  that 
can  be  found  explicitly.  For  Type  III  designs,  they  are: 


(i)  a  +  mb  =  -4p2  +  4p , 

(ii)  a  -  b  -  -2p  2  -  (4p-p 2 ) /m 


(3.2.3) 


It  is  easy  to  check  that  if  0  satisfies 

(1-2*) /2  <  0  <  {[l+m(m-l)/2]Js  -  l}/(m-l)  , 


(3.2.4) 


then  both  eigenvalues  are  between  -1  and  1,  so  that  the 
expansion  converges.  Expanding  (3.1.7)  gives 


8 

m 

s? 

V  -  *  *  ■ 

A  «  ■■ 

a  *  m-1{TV-T(2^N-p2N2)  (I-nT^T'T)  Y 


-  nr^T(2$N-£2N2)T'T(23N-£2N2)  (I-nT^T'T)Y  -  •••}.  (3.2.5) 

We  now  define  the  "modified  Papadakis  estimator"  as  the 
terms  of  zero  and  first  order  in  Q  -  2pN-p2N2  in  (3.2.5), 
namely 

oMP  «  m_1{TY-T(20N-02N2) (I-m"1T'T)Y)  .  (3.2.6) 

In  fact,  because  T(I-m  ^T'T)  =  0  we  can  simply  replace 
N2  by  N(2)  in  (3.2.6)  to  get 

SMP  -  m-1{TY-T (20N-02N ^2) )  (I-*m"’1T'T) Y)  .  (3.2.7) 

This  form  suggests  that  we  should  correct  the  mean  not 
only  by  the  effect  of  the  nearest  neighboring  plots  but 
also  by  the  effect  of  the  "lag  two"  neighboring  plots 
with  appropriate  coefficients,  as  shown  in  Table  1  . 


i-2  i-1  i  i+1  i+2 


Table  1.  Pattern  of  residual  correction  to  treat¬ 
ment  averages  for  the  one  dimention  "two 
sided"  error  structure. 
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This  estimator,  (3.2.7),  if  applied  iteratively  in  designs 
for  which  the  matrix  expansion  is  valid,  converges  for 
fixed  $  to  the  maximum  likelihood  estimator  (3.1.7). 

In  the  j iteration, 

aMPW  =m"1{TY-T(20N-02N(2))  (Y-T'aMP  (^_1)  )  }  (3.2.8) 

where  aMP^  is  m'^TY;  thus  aMP^  is  (3.2.7). 

3.3  An  Exact  Solution 

In  fact,  we  do  not  need  to  iterate  our  modified 
Papadakis  estimator  for  Type  III  designs  because  an  ex¬ 
plicit  solution  can  be  obtained  (below) ,  but  it  does  pro¬ 
vide  an  appropriate  parallel  to  the  ordinary  Papadakis 
estimator,  for  the  present  case.  For  Type  III  designs 
with  c  «  2,  TV~^T '  has  the  same  form  as  (2.3.8)  with 

a  =  m(l+2(32 )  ,  b  =  2(2g-02).  (3.3.1) 

Using  (2.4.2)  and  (2.4.3),  (3.1.5)  becomes 


V. V*  «\  .* ,  . 


r.' 


S  S 


X 


2 (m-1)  ( 2p-p 2 ) 

m (1+202 ) 


(1+20 


2)  1  - 


2(20-0') 

u+202) 


2  (2p-p) 

m(l+2027 


2g-02 _ 

m(l+202  )-2  (m-l)0  (2-0) 


x  TCI- (20N-02N2)  ]Yj.  (3.3.2) 

Note  that  TN^T'  is  a  t  by  t  matrix  whose  (s,w)t^1 
element  is  the  number  of  times  treatment  s  appears  adja¬ 
cent  but  one  to  treatment  w  in  the  design. 


3.4  Maximum  Likelihood  Estimator  for  p 

The  maximum  likelihood  equation  for  p  is  obtained 
in  several  stages.  We  first  find  the  maximum  likelihood 
estimator  for  o2, 


g2  =  (Y-T'S)  1  (I-0N)  2  (Y-T'ot) 


(3.4.1) 


then  we  differentiate  the  logarithm  of  Equation  (3.1.3) 
with  respect  to  p  and  set  the  result  equal  to  zero 
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to  give 

2) 

The  second  term  on  the  left  hand  side  of  (3.4.2)  has  an 
interesting,  continued  fraction  form.  First  of  all  the 
term  is  -2n  times  the  (i,i+l)  element  of  (I-f)N)”^.  This 
is  easy  to  observe  using  the  fact  that  the  derivative  of 
determinant  of  a  matrix  is  the  sum  of  the  determinants  of 
the  original  matrix  differentiated  column  by  column.  Here 


1  -p  0  •  •  •  0  0  -p 

-p  1  -p  •  •  •  o  0  0 

0-p  1  •  •  •  0  0  0 

3  ••••  ••• 

7p  ...  .  ... 

•  •  •  •  •  •  • 

0  0  0  •  •  •  1  -p  0 

0  0  0  •  •  •  -p  1  -p 

-p  0  0  •  •  •  0  -p  1 


n  (Y-T'a) '  (-2N+2gN2)  (Y-T'a)  ,  [ 3 1 1~pN|/9p]p=g  _  w  ... 

2  (Y-T'a) '  (I-20N+$ZNZ)  (Y-T'a)  |  I-pN| 


V  t.\ 


Si 


HI 

few 


a  V 

’•*  '/l 


i 


0  -1 


1  -p 


(3.4.3) 


1  -1 


Thus,  the  sum  of  these  determinants  is  -2n  times  the  co¬ 


factor  of  the  (i,i+l)  element  in  I-pN  and 


[3ll-pN|/3p3p^ 

J i-0n) 


(3.4.4) 


.  *  .  -1 

is  -2n  times  the  (i,i+l)  element  of  (I-pN)  .  Now 


abed 


babe 


c  b  a  b 


d  c  b  a 


(I  -  0N)-1= 


f  g  h 
e  f  g 


d  e 


e  d  c  b 


f  e  d  c 


g  f  e  d 
j  g  f  e 


(3.4.5) 


abed 


babe 


c  b  a  b 


d  c  b  a 


For  n  odd,  the  (i,i+l)  element  of  (I-£>N) 
which  we  illustrate  for  n  =  11: 


takes  a  form 


-p 


a2 


i  - 


*  2 


+  2a2  -  i 


i 


(3.4.6) 


1 


1-0 


In  general,  there  are  *i(n-3)  "p2  terms"  in  the  con¬ 
tinued  fraction.  For  n  even,  we  illustrate  for  n  =  12, 
for  which  the  required  expression  is 


-p 


a2 


+  202  -  1 


1  - 


1  - 


a2 


i  - 


l-2p 


it  2 


(3.4.7) 


In  general,  there  are  Mn-4)  ”£2  terms”  in  the  con¬ 
tinued  fraction,  not  counting  the  final  02  in  the  de¬ 
nominator  1-202  .  Thus,  an  explicit  solution  is  obtained 


by  substituting  the  a  in  Equation  (3. 3.4) into  Equation 
(3.4.2^  solving  for  £  ,  for  example,  using  the  Newton- 
Raphson  method,  and  then  obtaining  a  from  (3.3.4)  once 
0  is  found.  From  the  practical  point  of  view,  however, 
examination  of  the  likelihood  over  an  interval  of  p  , 
using  the  solution  for  a  in  terms  of  $  is  easier  than 
solving  explicitly  for  0  . 

(Note  that  if  all  $2  values  in  the  continued 
fractions  are  replaced  by  zeros,  we  obtain  0  ,  which 
thus  provides  an  approximate  solution.) 

3.5  Least  Squares  Solution  for  p 

The  least  squares  solution  for  p  is  the  3  that 
minimizes  the  weighted  residual  sum  of  squares.  This 
turns  out  to  be 


(Y-T 1 aLS) ,N(Y-T,gLS) 
(Y-T'aLS) ’N2 (Y-T'aLS) 


(3.5.1) 


which  is  similar  to  (2.3.5)  with  m  ^TY,  the  vector  of 
treatment  averages,  replaced  by  the  least  squares  esti¬ 
mator  of  a  .  The  latter  is  of  form  identical  to  a  , 


the  maximum  likelihood  solution,  but  with  the  least 
squares  0  inserted. 
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4.  TWO  DIMENSIONS,  "ONE-SIDED"  ERROR  STRUCTURE 

The  natural  extension  of  the  one  dimension,  first 
order,  one-sided  autoregressive  model  to  two  dimensions 
is  to  consider  the  model 


=  a  +  x, 


ID 


i  =  1, . 


j  =  1 


r 


s  =  l,...t,  (4.1) 


where 


i-l,j 


(4.2) 


and  e  ~  N(0,la2)  .  Let  cr  ■  n,  the  total  number  of 
observations.  The  likelihood  function  is 


(2Tro2)"n/2e3q>{-(2c2)’1  l  l  [x..  -  p(x.  .  .+x.  .  ,)]2)  = 

j— 1  i=i  AJ 


=  •  •  ■  =  (2ira2)_rv/2exp{-(2a2)''1[  (l+2p2)  £  [  x?  . 

j=l  i=l  1,3 


r  c 


r  c 


“  2P  I  I  xi  i  xi  i.i“2P  I  I  xi  -i  xi_i  -i 
j=l  1=1  X»D  X»D  J-  j=i  i-]^  ^D  1  !»D 

+  2p2  l  J  X.  .  .  X.  .  .  +  y  X?  +  f  X2  . 
j=l  i=l  1-1 '3  i,j-l  i,r  .£^  c,j 
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To  obtain  the  maximum  likelihood  estimator  for  a  we 

s 

first  replace  the  x..’s  in  (4.3)  by  Y.  .  -  a  .  The 

13  13  s 

appropriate  Jacobian  is  the  identity  matrix.  The  maxi 
mum  likelihood  estimator  is: 


«s  =  t"1{  E  yi-i  ”  l+lir  l  l  (yii '  «g> 

lCirj]=s  13  1+20  [i±l]*s  [j±l]=s  13  3 

+  T£>gT  ^  E  (yii““s} 

1+20  [i-l]=s  [j+l]=s  13  8 


+  1+202 


l  I  (yii-os)) 

Ci+l>s  [j-l]=s  13  s  J 


(4.4) 


where  l  is  the  number  of  times  treatment  s  was 

A 

applied,  i.e.,  the  maximum  likelihood  estimator  a  will 

s 

correct  the  average  of  all  the  yields  from  the  plots  re¬ 
ceiving  treatment  s  ,  by  the  residual  effects  of  the 
neighboring  plots  (i-l,j),  (i+1, j) ,  (i, j-1) ,  (i,j+l) 
with  coefficient  -  (l+2£2 )  ]  (plots  marked  with 

in  Table  2) ; 


and  by  those  of  plots  (i-l,j+l)  and  (i+l,j-l)  with 
coefficient  $2/[Jl  (1+202 )  ]  (plots  marked  with  x)  ;  how¬ 
ever  plots  (i-l,j-l)  and  (i+l,j+l)  are  ignored.  Equation 
(4.4)  can  be  rewritten  in  matrix  form  as 


a  *  l 


'Ht1 '  l’l'r(i«sTN ’  irSr  n<12>)  T‘l _1 
[z  -  (i«f  N ‘ n<12))]  y]  (4-5> 


where  T,  N,  Y,  a  defined  as  before  and  N 


is  the 


n  by  n  "semi  diagonal  neighboring"  matrix  symbolized  by 
x  in  Table  2. 


i-1  i  i+1 


j+1  x 


Table  2.  Pattern  of  residuals  corrections 
to  treatment  averages  for  the  two- 
dimens  i-.:  ns  "one-sided"  error 
structure . 


If  it  were  thought  reasonable  that  plots  ( i— 1 , j— 1 ) 
(i-l,j+l),  (i+l,j-l),  (i+l,j+l)  should  have  the  same 
influence  on  plot  ( i , j ) ,  this  "one-sided"  extension 
would  not  be  appropriate.  However,  one  blocking  factor 
could  be  time  or  influences  such  as  prevailing  wind 
direction  or  slope  of  the  ground  may  make  a  'one-sided' 
or  'semi-one-sided'  model  appropriate  in  some  contexts. 
Because  of  its  similarity  to  the  work  of  Section  5,  but 
with  a  redefinition  of  the  N'  matrix,  we  do  not 
deal  further  with  this  case  here. 
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5.  TWO  DIMENSIONS,  "TWO-SIDED”  ERROR  STRUCTURE 


5.1  The  Simultaneous  Autoregressive  Model 


Consider  a  lattice  consisting  of  a  finite  set  of 


sites,  each  site  having  associated  with  it  a  univariate 


random  variable.  To  quote  Besag  (1974,  p.  192):  "In 


most  ecological  applications,  the  sites  will  represent 


points  or  regions  in  the  Euclidean  plane  and  will  often 


be  subject  to  a  rigid  lattice  structure.  For  example. 


Cochran  (1936)  discussed  the  incidence  of  spotted  wilt 


over  a  rectangular  array  of  tomato  plants.  The  disease 


is  transmitted  by  insects  and,  after  an  initial  period 


of  time,  we  should  expect  to  observe  clusters  of  in¬ 


fected  plants . "  We  assume  the  model : 


y^  j  *  °s  +  xij'  is=lr2,...,c  j  =  l,2,...,r. 


s  —  l,2,...,t, 


(5.1.1) 


where 


xij  ■  p<xi-i,j  +  xi+i,j  +  xi,j-i  +  xi,j+i) +  eij-  <5-1-2) 


We  assume  that  e  ~  N(0,la2), 


S  V*. 


where  e  =  (en»e21'e31'  *'  *  'ecl'e12'  *  *  *'Ecr'  *  If  we  set 
n  -  rc,  the  likelihood  function  is 


r  c 


*  (2iro2)“n/2exp{-(2a2)"1C(l+4p2)C^1  ^  x2 

i=2  j*2  3 


c— 1  r-1 


o-l  r-1 


+  2p!iL  &*.&*** *&  jI2xi,J-lxi,3+l 


o-l  r-1 


c— 1  r-1 


+  4p2  L-2  jLXi.jXi+1.3+1+4p!  il2 


c-1  r-1 


c— 1  r-1 


4f>  1=2  jL*1**1*1' j  "  4P  iL  jI2Xi^Xi'j+l]}  * 


(5.1.3) 


The  maximum  likelihood  estimator  of  a  is 

s 


“s  *  1_1{  £  yij'I^iF  *  J  (yij'sij’ 

s  Hi,  j  ]=s  13  1  *p  Ci±l]=s  [j±l]=s  J  J 


+  -ilfp-  £  ^  (yii-sn' 

1+4p  [i+i]=s  [j±i]-s  13  13 


+  l+40s 


C i± 2 ]=s  [ j± 2 ]=s 


I  (y .  .-a .  .  A 

i±  2 l  =  s  13  13  i 


(5.1.4) 


where  l  is  the  number  of  times  treatment  s  was 


applied.  Here,  the  average  of  all  the  plots  receiving 
treatment  s  will  be  corrected  for  the  residual  effects 
of  the  neighboring  plots  (i-l,j),  (i+l,j),  (i,j-l),  (i,j+l) 
with  coefficient  -  20/C l (1+402 ) ]  (shown  as  •  in  Table 
3);  for  those  of  plots  (i-1, j-1) ,  (i-l,j+l),  (i+l,j-l), 
(i+1, j+1)  with  coefficient  202/Ci (1+402 ) ]  (shown  as  x) ; 
and  for  those  of  plots  (i-2,j),  (i+2,j),  (i,j-2),  (i,j+2) 
with  coefficient  02/C i ( 1+40 2 )  ]  (shown  as  *)  . 


Table  3.  Pattern  of  residuals  corrections  to 
treatment  averages  for  the  two- 
dimensional,  !' two-sided"  error 
structure. 


We  now  recast  Equation  (5.1.4)  into  matrix  form.  The 
likelihood  function  and  the  maximum  likelihood  estimator 
for  a  are  exactly  the  forms  of  Equations  (3.1.3)  and 
(3.1.5),  but  with  redefined  N  and  T.  For  example, 
for  the  5x5  knight's  move  Latin  square  given  by 


(5.1.5) 


we  have 


T 

5x25 


1000000100000010100000010 
0100000010100000010000001 
0010000001010000001010000  (5.1.6) 

0001010000001000000101000 
0000101000000101000000100 


The  n  by  n  neighbor-specification  matrix  N  is  now 
different  from  that  of  Chapters  2  and  3.  Here 


N  =  Nc  ®  Ir 
25x25  5x5  5x3 


[c  ®  N* 


(5.1.7) 


,*>v>  Wm 


>  VV  VW.v,  /.  < ,  *.  •.  -  .  ■  .  ■  .  «  .  •  ■.  s  y.  .'.A  .  •  . 


.  •*.  •*  •  • 


where 


0  10  0  1 
10  10  0 
0  10  10 
0  0  10  1 
10  0  10 

and  ®  is  the  Kronecker  product.  In  general,  Nr  and 
* 

N  are  the  r  by  r  and  c  by  c  neighbor-specifica- 

C 

tion  matrices  for  rows  and  columns  respectively.  We  can 
write 

N2  =  N(2)  +  2N(12)  +  41  (5.1.9) 

where,  with  a  slight  change  of  notation  from  the  one- 

(12) 

dimensional  case,  which  offers  no  "diagonal  plots",  N 

is  the  "diagonal  neighbor-specification"  matrix,  symbol- 

(2) 

ized  by  x  in  Table  3,  and  N  is  the  "lag  two" 
neighbor-specification  matrix,  symbolized  by  *  in  Table 
3.  The  •  symbolizes  N  .  Then, 

V”1  -  (I-0N)  2  ■  (1+402)!  -  23N  +  2£2N(12)  +  £2N(2).  (5.1.10) 


We  thus  obtain,  from  (3.1.5)  and  (5.1.10), 


rT  (  28  2i!\  ,,02)  S! 

[I_  \l+49!  N"l+4$!  N  1+49' 


n(2,)J  y|  ,  (5 


.l.ii) 


which  is  identical  to  (5.1.4). 


5.2  The  Modified  Papadakis  Estimator 

The  maximum  likelihood  estimator  can  be  expanded  if 
and  only  if  all  the  eigenvalues  of  8.  *T  (2f)N-(52N2 ) T'  are 
smaller  than  one  in  absolute  value.  For  example,  for  two 
dimensional  designs  derived  by  making  a  row  by  row  cyclic 
permutation  mt  times  of  the  Type  III  designs  of  Williams 
(1952,  p.  159),  £-1T(23N-02N2)T'  is  a  pattern  matrix  of 
the  form  of  (2.3.8)  with 


a  =  -8pz  , 
b  =  8 ( p— p  2 ) /m  , 


(5.2.1) 


with  the  two  distinct  eigenvalues 
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(i) 

a  +  mb  =  -8p2 

+  8 (p-p 2  )  , 

(ii) 

a  -  b  =  -8p2 

-  8 (p-p2 )/m  . 

If  0  satisfies 

(1-2*) /4  <  $  <  [ { 1-Hn (m-1) /2) **-1  ]/{ 2 (m-1) } ,  (5.2.3) 

the  expansion  is  valid  and  takes  the  form 

a  =  *_1{TY-T(20N-02N2)  (IhTV'DY 

-  i”^2£N-02N2)T,T(20N-02N2)  (I-JTVtJY- •••  }.  (5.2.4) 

Deletion  of  all  terms  in  0  whose  powers  exceed  one  gives 
an  estimator  of  the  Papadakis  type  for  the  two  dimensional 
case 


£“1{TY  -  20TN(I-fc_1T'T) Y)  (5.2.5) 

which  may  be  compared  with  (2.3.4).  If  we  retain  terms 
of  order  zero  and  one  in  Q  =  2f5N-(52N2,  so  that  cor¬ 
rections  for  diagonal  and  "lag  two"  plots  are  included, 
we  can  define  a  "modified  Papadakis  estimator”: 


1  {TY  -  T  (20N  -  j$zN2  )  (I  -  l  T'T)  Y> 


(5.2.6) 


and  we  achieve  the  same  phenomena  as  in  the  one  dimen¬ 
sional  case  of  reproducing  the  term  Jt"^T(2£N-02N2  )T' 
again  and  again.  If  Bartlett's  (1978)  suggestion  of 
iterating  is  followed  for  the  "modified  Papadakis  esti¬ 
mator"  (5.2.6) ,  convergence  to  the  maximum  likeli¬ 
hood  estimator  occurs,  for  0  fixed,  for  designs  for 
which  the  matrix  expansion  is  valid.  Because  T(I-m  ^T'T) 
0,  we  can  also  write  this  "modified  Papadakis  estimator" 
in  the  form 


(5.2.7) 


i-1{TY-T (2$N-20*N <12)  -02N (2) )  (I-t~^T'T) Y} 


which  directly  reflects  the  pattern  of  Table  3. 


5.3  An  Exact  Solution 


Iteration  is  not  actually  needed  for  designs  of  the 
type  described  above,  however.  As  in  the  one-dimensional 
case,  we  can  again  get  an  explicit  solution  for  cyclically 
permuted  (t  times,  or  mt  times,  with  appropriate 
bordering  treatments,  different  for  the  two  cases)  Type 
III  designs  without  resorting  to  the  Papadakis  estimator, 


.’.y ’A  -■  «*.  .-.y  .v .■.-y.-  .vy  v  .'.v.v.v.v,'.' 

*  *  *  i*»  *  »**#’.  •  ■  *  *  *  h  *  •  *  ,  *  j  *  »  •  •  *  •  *  B  ••  *  *  i>  *  »  »ji  ■  „  ‘kN  ,  wV.%  V  * 


by  using  pattern  results  to  invert  the  matrix  TV  T' 
that  has  the  form  of  (2.3.8)  with 

a  =  m2  (m+1)  (1+802)  » 

b  =  -m(m+l)  (80-8g2)  .  (5.3.1) 

For  example,  the  solution  for  Type  III  designs  cyclically 
permuted  mt  times  is 


1  (m~l) (8p-8p2 ) 

1 _ m(l+8p2) 


Note  that  TN^^T'  is  a  t  by  t  matrix  whose  (s,w)th 

element  is  the  number  of  times  treatment  s  appears  on 

(2) 

the  diagonal  of  treatment  w  in  the  design,  and  TN  T* 
is  a  t  by  t  matrix  whose  (s,w)  element  is  the  number  of 
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times  treatment  s  appears  adjacent  but  one  to  treatment 
w  in  the  design. 

5.4  Maximum  Likelihood  Estimator  for  p 

The  maximum  likelihood  equation  for  £  is  obtained 
as  in  section  (3.4)(pp.  32-35)  with  symbols  appropriately 
redefined  for  two  dimensions.  We  now  need  -4n  times  the 
(i,i+l)  element  of  B-1  =  (I-0N)"1.  Because  of  the 
form  of  B,  B-1  will  have  a  special  form.  The  first 
column  of  B-1,  b^,  say,  consists  of  a  small  subset 
containing  q,  say,  elements  which  recur  throughout  B 
The  n  equations  B^  =  (1,0,0, ...  ,0)  '  reduce  to  q 
equations  B*b*  «  (1,0,0,. ..,0)'  where  B*  is  qxq  and 
both  b*  and  the  right-hand-side  are  qxl  vectors.  The 
element  we  want  is  then  (B* ) 7^  where  j  is  the  position 
of  the  required  element  in  b*,  this  position  being  our 
arbitrary  choice.  The  value  of  q  is 

RU 

q  =  \  i  ,  (5.4.1) 

i=l 

where 

RU  =  RU{»s  (n*5  +  1)  }  , 


(5.4.2) 


and  RU (0 )  means  "round  up  0  to  the  next  integer  if 
it  is  not  already  an  integer" .  A  specific  example  for 
such  a  procedure  is  given  in  Appendix  2.  As  before,  we 
now  solve  iteratively  for  p  in  Equation  (5.1.11),  or 
(5.3.2) . 

The  least  squares  solution  for  p  is  similar  to 
(3.5.1)  with  l  replacing  m,  and  with  N  appropriately 


redefined. 


6.  SIMULATION  STUDY 


A  simulation  study  to  compare  the  different  estima¬ 
tors  described  in  Chapters  2,  3,  and  5  was  conducted  both 
for  one  and  two  dimensions. 

6.1  One  Dimension 

In  one  dimension,  four  estimation  procedures  were 
compared : 

(i)  The  maximum  likelihood  estimator  for  the 
simultaneous  autoregressive  model  (Section 

3.1)  . 

(ii)  The  modified  Papadakis  estimator  (Section 

3.2)  . 

(iii)  The  maximum  likelihood  estimator  for  the 
first  order  autoregressive  model  (Section 
2.1) . 

(iv)  The  Papadakis  estimator  (section  2.2). 

The  design  assumed  was  the  Type  III  design 

(1, 2, 3, 4) (2, 1,3, 4) (1,3,2, 4) .  (6.1.1) 

The  observations  were  generated  as  follows: 


(1)  The  model  assumed  was 

Y  *  T'o  +  X/  (6.1.2) 

where 

X  -  pNX  +  e  , 

so  e  ■  (I  -  pN) x  =—*>  X  *  (I  -  pN) -1e  . 

(2)  Using  generated  normal  random  variables  with 
zero  mean  and  three  different  standard  deviations 
(0.01,0.1,1)  X  was  found  for  fixed  p  . 

(3)  On  the  X  was  superimposed  a  treatment  effect 
as  required  by  the  design  (6.1.1)  to  get  the 
observations  Y  . 

(4)  The  above  scheme  was  used  to  construct  twenty 
data  sets  for  various  combinations  of  p 


The  four  different  estimation  procedures  were  then 
applied  to  all  20  data  sets. 

In  Table  4,  the  maximum  likelihood  estimator  for 
p  is  given  as  well  as  the  weighted  residuals  sum  of 
squares  for  the  different  models  and  methods.  Note 
that,  for  the  modified  Papadakis  estimator,  p  was  de¬ 
rived  from  the  likelihood  function  of  the  simultaneous 
autoregressive  model  with  the  modified  Papadakis  esti¬ 
mator  for  a  inserted;  this  explains  the  difference 
in  the  p  for  the  maximum  likelihood  estimator  and  the 
p  for  the  modified  Papadakis  estimator.  The  weighted 
residuals  sum  of  squares  was  calculated  as 

(Y-T'SWY- T'a)  (6.1.3) 

where 

(1)  V  ^  *  (I  -  pN) 2  for  the  estimator  from 

the  simultaneous  autore¬ 
gressive  model. 

(2)  V  ^  =  (I  -  pN) 2  for  the  modified  Papadakis 

estimator. 


(3)  V  *  I  -  YN  for  the  estimation  from  the 

first  order  autoregressive 
model. 

(4)  V  1  =  I  -  JjfjN  for  the  Papadakis  estimator. 

We  observe,  from  Table  4,  that  overall  the  simul 
taneous  autoregressive  model  does  better,  as  might  be 
expected.  Only  for  certain  p  values,  as  tabulated, 
do  the  other  methods  provide  close  estimates  of  p  . 
The  weighted  residual  sum  of  squares  is  smallest  for 
the  simultaneous  autoregressive  model  everywhere  ex¬ 
cept  for  simulation  7  ,  and  19  where  the  modified 
Papadakis  estimator  is  slightly  better. 


Table  4.  (Continued) 


6 . 2  Two  Dimensions 

For  two  dimensions,  three  estimators  were  compared 

(i)  The  maximum  likelihood  estimator  from  the 
simultaneous  autoregressive  model  (Section 

5.1) . 

(ii)  The  modified  Papadakis  estimator  (Section 

5.2)  . 

(iii)  The  extension  of  the  Papadakis  estimator 
into  two  dimensions  (Section  5.2). 

The  design  used  was  (6.1.1)  row  by  row  cyclically 
permuted  12  times.  Thirteen  different  data  sets  were 
generated  in  a  manner  similar  to  that  described  for  the 
one  dimensional  case.  Table  5  is  constructed  in  the 
same  manner  as  was  Table  4  for  the  one  dimensional  case, 
with  the  obvious  extension  into  two  dimensions. 

Table  5  shows  the  same  characteristics  for  two 
dimensions  that  Table  4  showed  for  one  dimension  ex¬ 
cept  that  the  comparison  flatters  the  simultaneous  auto¬ 
regressive  model  even  more.  Note  that  the  data  in 
Table  5  involves  144  observations,  compared  with  the 


12  used  in  Table  4 


7 .  REMARK 


Expansion  of  the  MLE  equations  for  any  given  model 
set-up  may  or  may  not  be  valid.  If  it  is,  the  "zero  and 
first  order"  terms  provide  an  appropriate  Papadakis-type 
estimator  for  the  problem  studied,  and  iteration  on  this 
estimator  converges  to  the  MLE.  Thus  the  point  made  by 
Wilkinson  et  al.  (1983)  that  the  iterated  Papadakis 
estimator  will  provide  a  positively  biased  treatment 
F-ratio  is  puzzling.  A  possible  explanation  is  that  the 
Papadakis-type  estimators  used  by  Wilkinson  et  al.  (1983) 
may  not  be  the  appropriate  Papadakis-type  estimator  pro¬ 
vided  by  the  first  term  in  a  valid  matrix  expansion  under 
the  model  assumed.  Alternatively,  such  an  expansion  may 
not  be  valid.  These  issues  appear  to  need  further 
exploration. 


Appendix  1.  Some  Useful  One  and  Two  Dimensional 
Neighbor  Balanced  Designs  Reproduced 
from  Williams  1952  and  Freeman  1979. 
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(See  Section  2.3) 

One  Dimension:  Type  11(a)  Designs  with  c  =  2 

(Williams  1952). 

In  Type  11(a)  designs,  each  treatment  occurs  n 
times  and  equally  often,  c  times,  say,  next  to  every 
other  treatment. 

m  -  4  (1,2, 3, 4, 5)  (3, 4, 1,5, 2)  (4, 5, 3, 1,2)  (3, 1,4, 2, 5) 

m  =  5  (1,2, 3, 4, 5, 6)  (2, 5, 3, 6, 4,1)  (5, 3, 1,6, 4, 2)  (5, 1,6, 3, 2, 4) 

(3. 1.4. 5. 6. 2) 

6  (1,2, 3, 4, 5, 6, 7) (5, 3, 6, 4, 7, 2,1)  (4, 2, 5, 1,6, 3, 7) 

(1.3. 5. 4. 7. 6. 2)  (5, 6, 4, 3, 2, 7,1)  (6, 2, 4, 1,5, 7, 3) 

m  ■  7  (1, 2, 3, 4, 5.,6, 7, 8)  (6, 4, 2, 5, 1,3, 8, 7)  (5, 3, 6, 2, 7, 4, 8,1) 

(7. 3. 4. 1.6. 5. 8. 2)  (6, 3, 7, 4, 8, 1,2, 5) 

(3, 2, 4, 5, 8, 6, 7,1)  (6, 4, 1,3, 8, 2, 7, 5) 


(Williams  1952). 

Type  III  designs  satisfy  the  condition  of  Type  11(a) 
designs;  also  each  treatment  occurs  equally  often  adja¬ 
cent  but  one  to  every  other  treatment. 

ra  »  2  (1,2,3) (1,2,3) 

m=3  (1,2, 3, 4,) (2, 1,3, 4) (1,3, 2, 4) 

m«  4  (1,2, 3, 4, 5) (2, 4, 1,5, 3) (1,4, 5, 3, 2) (5, 1,2, 4, 3) 

Two  Dimensions:  Complete  Latin  Squares 

(Freeman  1979). 

Complete  Latin  squares  are  designed  such  that  any 
two  treatments  occur  next  to  each  other  once  in  a  row 
and  once  in  a  column.  (Note  that  the  designs  are  not 
assumed  to  be  torus  designs.) 

12  3  4 

3  14  2 

m  =  4 

2  4  13 


4  3  2  1 


5? 


Appendix  2.  An  Example  of  the  Derivation  of  the 
Maximum  Likelihood  Estimator  for  p  . 
(See  Section  5.4) 


Consider  the  5  by  5  Latin  square: 

1  2  3  4  5 

2  3  4  5  1 

3  4  5  1  2 

4  5  12  3 

5  12  3  4 


B*1*  (I-pN)"1 
25x25 


R 

1 

s  1 

T 

1 

T 

1 

s  1 

—  - 

-1 

-  -1 

-  - 

—  — 

-  — 

S 

1 

R  i 

S 

1 

T 

1 

T 

—  - 

— 

— ■  - 

—  — 

-  - 

T 

1 

s  , 

R 

1 

S 

1 

T 

T 

1 

T  1 

S 

1 

R 

I 

S 

—  - 

-1 

- 1 

—  - 

—  — 

-  - 

S 

1 

T  i 

T 

t 

S 

1 

R 
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The  25  equations  Bb^  =  (1,0,0,...,0)  now  reduce  to 


q  s  6  equations  B*b*  —  (1,0,0»0,0,0)  where 

6*6  6*6 


A  A 

-p  -p  0 


-p  1-p  0  -2p  0 


1  -2p  0 


A  A  ,  A  A 

-p  -p  1-p  -p 


0  -2p  1-20 


and  b*  =  (a,b,c,d,e,f )  ' .  We  need  only  to  solve  for 


the  portion  of  the  solution  involving  element  b  of 


b*  . 


1 
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